home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / ranlib / tstmid.for < prev    next >
Text File  |  1996-07-19  |  13KB  |  449 lines
  1.       SUBROUTINE stat(x,n,av,var,xmin,xmax)
  2. C**********************************************************************
  3. C
  4. C     SUBROUTINE STAT( X, N, AV, VAR)
  5. C
  6. C               compute STATistics
  7. C
  8. C
  9. C                              Function
  10. C
  11. C
  12. C     Computes AVerage and VARiance of array X(N).
  13. C
  14. C**********************************************************************
  15. C     .. Scalar Arguments ..
  16.       REAL av,var,xmax,xmin
  17.       INTEGER n
  18. C     ..
  19. C     .. Array Arguments ..
  20.       REAL x(n)
  21. C     ..
  22. C     .. Local Scalars ..
  23.       REAL sum
  24.       INTEGER i
  25. C     ..
  26. C     .. Intrinsic Functions ..
  27.       INTRINSIC real
  28. C     ..
  29. C     .. Executable Statements ..
  30.       xmin = x(1)
  31.       xmax = x(1)
  32.       sum = 0.0
  33.       DO 10,i = 1,n
  34.           sum = sum + x(i)
  35.           IF (x(i).LT.xmin) xmin = x(i)
  36.           IF (x(i).GT.xmax) xmax = x(i)
  37.    10 CONTINUE
  38.       av = sum/real(n)
  39.       sum = 0.0
  40.       DO 20,i = 1,n
  41.           sum = sum + (x(i)-av)**2
  42.    20 CONTINUE
  43.       var = sum/real(n-1)
  44.       RETURN
  45.  
  46.       END
  47.       PROGRAM tstall
  48.       IMPLICIT LOGICAL (q)
  49. C     Interactive test for PHRTSD
  50. C     .. Parameters ..
  51.       INTEGER mxwh
  52.       PARAMETER (mxwh=10)
  53. C     ..
  54. C     .. Local Scalars ..
  55.       REAL av,avtr,var,vartr,xmin,xmax
  56.       INTEGER i,is1,is2,itmp,iwhich,j,mxint,nperm,nrep,ntot
  57.       CHARACTER type*4,phrase*100
  58. C     ..
  59. C     .. Local Arrays ..
  60.       REAL array(1000),param(3)
  61.       INTEGER iarray(1000),perm(500)
  62. C     ..
  63. C     .. External Functions ..
  64.       REAL genbet,genchi,genf,gennch,gennf,genunf
  65.       INTEGER ignuin,myhand
  66.       EXTERNAL genbet,genchi,genf,gennch,gennf,genunf,ignuin,myhand
  67. C     ..
  68. C     .. External Subroutines ..
  69.       EXTERNAL genprm,phrtsd,setall,stat,trstat
  70. C     ..
  71. C     .. Equivalences ..
  72.       EQUIVALENCE (array,iarray)
  73. C     ..
  74. C     .. Executable Statements ..
  75.       WRITE (*,9000)
  76.  
  77.  9000 FORMAT (' Tests most generators of specific distributions.'/
  78.      +       ' Generates 1000 deviates: reports mean and variance.'/
  79.      +       ' Also reports theoretical mean and variance.'/
  80.      +       ' If theoretical mean or var doesn''t exist prints -1.'/
  81.      +       ' For permutations, generates one permutation of 1..n'/
  82.      +       '     and prints it.'/
  83.      +       ' For uniform integers asks for upper bound, number of'/
  84.      +       '     replicates per integer in 1..upper bound.'/
  85.      +       '     Prints table of num times each integer generated.')
  86. C
  87. C     Menu for choosing tests
  88. C
  89.    10 WRITE (*,9010)
  90.  
  91.  9010 FORMAT (' Enter number corresponding to choice:'/
  92.      +       '      (0) Exit this program'/
  93.      +       '      (1) Generate Chi-Square deviates'/
  94.      +       '      (2) Generate noncentral Chi-Square deviates'/
  95.      +       '      (3) Generate F deviates'/
  96.      +       '      (4) Generate noncentral F  deviates'/
  97.      +       '      (5) Generate random permutation'/
  98.      +       '      (6) Generate uniform integers'/
  99.      +       '      (7) Generate uniform reals'/
  100.      +       '      (8) Generate beta deviates'/
  101.      +       '      (9) Generate binomial outcomes'/
  102.      +       '     (10) Generate Poisson ourcomes'/)
  103.  
  104.       READ (*,*) iwhich
  105.       IF (.NOT. (iwhich.LT.0.OR.iwhich.GT.mxwh)) GO TO 20
  106.       WRITE (*,*) ' Choices are 1..',mxwh,' - try again.'
  107.       GO TO 10
  108.  
  109.    20 IF (iwhich.EQ.0) STOP ' Normal termination rn tests'
  110.       WRITE (*,*) ' Enter phrase to initialize rn generator'
  111.       READ (*,'(a)') phrase
  112.       CALL phrtsd(phrase,is1,is2)
  113.       CALL setall(is1,is2)
  114.  
  115.       IF ((1).NE. (iwhich)) GO TO 40
  116. C
  117. C     Chi-square deviates
  118. C
  119.       type = 'chis'
  120.       WRITE (*,*) ' Enter (real) df for the chi-square generation'
  121.       READ (*,*) param(1)
  122.       DO 30,i = 1,1000
  123.           array(i) = genchi(param(1))
  124.    30 CONTINUE
  125.       CALL stat(array,1000,av,var,xmin,xmax)
  126.       CALL trstat(type,param,avtr,vartr)
  127.       WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
  128.  
  129.  9020 FORMAT (' Mean Generated: ',T30,G15.7,5X,'True:',T60,
  130.      +       G15.7/' Variance Generated:',T30,G15.7,5X,'True:',T60,
  131.      +       G15.7/' Minimun: ',T30,G15.7,5X,'Maximum:',T60,G15.7)
  132.  
  133.       GO TO 280
  134.  
  135.    40 IF ((2).NE. (iwhich)) GO TO 60
  136.  
  137. C
  138. C     Noncentral Chi-square deviates
  139. C
  140.       type = 'ncch'
  141.       WRITE (*,*) ' Enter (real) df'
  142.       WRITE (*,*) '       (real) noncentrality parameter'
  143.       READ (*,*) param(1),param(2)
  144.       DO 50,i = 1,1000
  145.           array(i) = gennch(param(1),param(2))
  146.    50 CONTINUE
  147.       CALL stat(array,1000,av,var,xmin,xmax)
  148.       CALL trstat(type,param,avtr,vartr)
  149.       WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
  150.       GO TO 280
  151.  
  152.    60 IF ((3).NE. (iwhich)) GO TO 80
  153.  
  154. C
  155. C     F deviates
  156. C
  157.       type = 'f'
  158.       WRITE (*,*) ' Enter (real) df of the numerator'
  159.       WRITE (*,*) '       (real) df of the denominator'
  160.       READ (*,*) param(1),param(2)
  161.       DO 70,i = 1,1000
  162.           array(i) = genf(param(1),param(2))
  163.    70 CONTINUE
  164.       CALL stat(array,1000,av,var,xmin,xmax)
  165.       CALL trstat(type,param,avtr,vartr)
  166.       WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
  167.       GO TO 280
  168.  
  169.    80 IF ((4).NE. (iwhich)) GO TO 100
  170.  
  171.  
  172. C
  173. C     Noncentral F deviates
  174. C
  175.       type = 'ncf'
  176.       WRITE (*,*) ' Enter (real) df of the numerator'
  177.       WRITE (*,*) '       (real) df of the denominator'
  178.       WRITE (*,*) '       (real) noncentrality parameter'
  179.       READ (*,*) param(1),param(2),param(3)
  180.       DO 90,i = 1,1000
  181.           array(i) = gennf(param(1),param(2),param(3))
  182.    90 CONTINUE
  183.       CALL stat(array,1000,av,var,xmin,xmax)
  184.       CALL trstat(type,param,avtr,vartr)
  185.       WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
  186.       GO TO 280
  187.  
  188.   100 IF ((5).NE. (iwhich)) GO TO 140
  189.  
  190. C
  191. C     Random permutation
  192. C
  193.   110 WRITE (*,*) ' Enter size of permutation'
  194.       READ (*,*) nperm
  195.       IF (.NOT. (nperm.LT.1.OR.nperm.GT.500)) GO TO 120
  196.       WRITE (*,*) ' Permutation size must be between 1 and 500 ',
  197.      +  '- try again!'
  198.       GO TO 110
  199.  
  200.   120 WRITE (*,*) '       Random Permutation Generated - Size',nperm
  201.       DO 130,i = 1,500
  202.           perm(i) = i
  203.   130 CONTINUE
  204.       CALL genprm(perm,nperm)
  205.       WRITE (*,*) ' Perm Generated'
  206.       WRITE (*,'(20I4)') (perm(i),i=1,nperm)
  207.       GO TO 280
  208.  
  209.   140 IF ((6).NE. (iwhich)) GO TO 170
  210.  
  211. C
  212. C     Uniform integer
  213. C
  214.       WRITE (*,*) ' Enter maximum uniform integer'
  215.       READ (*,*) mxint
  216.       WRITE (*,*) ' Enter number of replications per integer'
  217.       READ (*,*) nrep
  218.       DO 150,i = 1,1000
  219.           iarray(i) = 0
  220.   150 CONTINUE
  221.       ntot = mxint*nrep
  222.       DO 160,i = 1,ntot
  223.           itmp = ignuin(1,mxint)
  224.           iarray(itmp) = iarray(itmp) + 1
  225.   160 CONTINUE
  226.       WRITE (*,*) '         Counts of Integers Generated'
  227.       WRITE (*,'(20I4)') (iarray(j),j=1,mxint)
  228.       GO TO 280
  229.  
  230.   170 IF ((7).NE. (iwhich)) GO TO 190
  231.  
  232. C
  233. C     Uniform real
  234. C
  235.       type = 'unif'
  236.       WRITE (*,*) ' Enter Low then High bound for uniforms'
  237.       READ (*,*) param(1),param(2)
  238.       DO 180,i = 1,1000
  239.           array(i) = genunf(param(1),param(2))
  240.   180 CONTINUE
  241.       CALL stat(array,1000,av,var,xmin,xmax)
  242.       CALL trstat(type,param,avtr,vartr)
  243.       WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
  244.       GO TO 280
  245.  
  246.   190 IF ((8).NE. (iwhich)) GO TO 210
  247.  
  248. C
  249. C     Beta deviate
  250. C
  251.       type = 'beta'
  252.       WRITE (*,*) ' Enter A, B for Beta deviate'
  253.       READ (*,*) param(1),param(2)
  254.       DO 200,i = 1,1000
  255.           array(i) = genbet(param(1),param(2))
  256.   200 CONTINUE
  257.       CALL stat(array,1000,av,var,xmin,xmax)
  258.       CALL trstat(type,param,avtr,vartr)
  259.       WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
  260.       GO TO 280
  261.  
  262.   210 IF ((9).NE. (iwhich)) GO TO 240
  263. C
  264. C     Binomial outcomes
  265. C
  266.       type = 'bin'
  267.       WRITE (*,*) ' Enter number of trials, Prob event for ',
  268.      +  'binomial outcomes'
  269.       READ (*,*) ntry,pevt
  270.       DO 220,i = 1,1000
  271.           iarray(i) = ignbin(ntry,pevt)
  272.   220 CONTINUE
  273.       DO 230,i = 1,1000
  274.           array(i) = iarray(i)
  275.   230 CONTINUE
  276.       CALL stat(array,1000,av,var,xmin,xmax)
  277.       param(1) = ntry
  278.       param(2) = pevt
  279.       CALL trstat(type,param,avtr,vartr)
  280.       WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
  281.       GO TO 280
  282.  
  283.   240 IF ((10).NE. (iwhich)) GO TO 270
  284. C
  285. C     Poisson outcomes
  286. C
  287.       type = 'pois'
  288.       WRITE (*,*) ' Enter mean for Poisson generation'
  289.       READ (*,*) param(1)
  290.       DO 250,i = 1,1000
  291.           iarray(i) = ignpoi(param(1))
  292.   250 CONTINUE
  293.       DO 260,i = 1,1000
  294.           array(i) = iarray(i)
  295.   260 CONTINUE
  296.       CALL stat(array,1000,av,var,xmin,xmax)
  297.       CALL trstat(type,param,avtr,vartr)
  298.       WRITE (*,9020) av,avtr,var,vartr,xmin,xmax
  299.       GO TO 280
  300.  
  301.   270 CONTINUE
  302.   280 GO TO 10
  303.  
  304.       END
  305.       SUBROUTINE trstat(type,parin,av,var)
  306.       IMPLICIT INTEGER (i-n),REAL (a-h,o-p,r-z),LOGICAL (q)
  307. C**********************************************************************
  308. C
  309. C     SUBROUTINE TRSTAT( TYPE, PARIN, AV, VAR )
  310. C               TRue STATistics
  311. C
  312. C     Returns mean and variance for a number of statistical distribution
  313. C     as a function of their parameters.
  314. C
  315. C
  316. C                              Arguments
  317. C
  318. C
  319. C     TYPE --> Character string indicating type of distribution
  320. C             'chis' chisquare
  321. C             'ncch' noncentral chisquare
  322. C             'f'    F (variance ratio)
  323. C             'ncf'  noncentral f
  324. C             'unif' uniform
  325. C             'beta' beta distribution
  326. C                         CHARACTER*(4) TYPE
  327. C
  328. C     PARIN --> Array containing parameters of distribution
  329. C              chisquare
  330. C               PARIN(1) is df
  331. C              noncentral chisquare
  332. C               PARIN(1) is df
  333. C               PARIN(2) is noncentrality parameter
  334. C              F (variance ratio)
  335. C               PARIN(1) is df numerator
  336. C               PARIN(2) is df denominator
  337. C              noncentral F
  338. C               PARIN(1) is df numerator
  339. C               PARIN(2) is df denominator
  340. C               PARIN(3) is noncentrality parameter
  341. C              uniform
  342. C               PARIN(1) is LOW bound
  343. C               PARIN(2) is HIGH bound
  344. C              beta
  345. C               PARIN(1) is A
  346. C               PARIN(2) is B
  347. C                         REAL PARIN(*)
  348. C              binonial
  349. C               PARIN(1) is Number of trials
  350. C               PARIN(2) is Prob Event at Each Trial
  351. C              poisson
  352. C               PARIN(1) is Mean
  353. C
  354. C     AV <-- Mean of specified distribution with specified parameters
  355. C                         REAL AV
  356. C
  357. C     VAR <-- Variance of specified distribution with specified paramete
  358. C                         REAL VAR
  359. C
  360. C
  361. C                              Note
  362. C
  363. C
  364. C     AV and Var will be returned -1 if mean or variance is infinite
  365. C
  366. C**********************************************************************
  367. C     .. Scalar Arguments ..
  368.       REAL av,var
  369.       CHARACTER type* (4)
  370. C     ..
  371. C     .. Array Arguments ..
  372.       REAL parin(*)
  373. C     ..
  374. C     .. Local Scalars ..
  375.       REAL a,b,range
  376. C     ..
  377. C     .. Executable Statements ..
  378.       IF (('chis').NE. (type)) GO TO 10
  379.       av = parin(1)
  380.       var = 2.0*parin(1)
  381.       GO TO 170
  382.  
  383.    10 IF (('ncch').NE. (type)) GO TO 20
  384.       a = parin(1) + parin(2)
  385.       b = parin(2)/a
  386.       av = a
  387.       var = 2.0*a* (1.0+b)
  388.       GO TO 170
  389.  
  390.    20 IF (('f').NE. (type)) GO TO 70
  391.       IF (.NOT. (parin(2).LE.2.0001)) GO TO 30
  392.       av = -1.0
  393.       GO TO 40
  394.  
  395.    30 av = parin(2)/ (parin(2)-2.0)
  396.    40 IF (.NOT. (parin(2).LE.4.0001)) GO TO 50
  397.       var = -1.0
  398.       GO TO 60
  399.  
  400.    50 var = (2.0*parin(2)**2* (parin(1)+parin(2)-2.0))/
  401.      +      (parin(1)* (parin(2)-2.0)**2* (parin(2)-4.0))
  402.    60 GO TO 170
  403.  
  404.    70 IF (('ncf').NE. (type)) GO TO 120
  405.       IF (.NOT. (parin(2).LE.2.0001)) GO TO 80
  406.       av = -1.0
  407.       GO TO 90
  408.  
  409.    80 av = (parin(2)* (parin(1)+parin(3)))/ ((parin(2)-2.0)*parin(1))
  410.    90 IF (.NOT. (parin(2).LE.4.0001)) GO TO 100
  411.       var = -1.0
  412.       GO TO 110
  413.  
  414.   100 a = (parin(1)+parin(3))**2 + (parin(1)+2.0*parin(3))*
  415.      +    (parin(2)-2.0)
  416.       b = (parin(2)-2.0)**2* (parin(2)-4.0)
  417.       var = 2.0* (parin(2)/parin(1))**2* (a/b)
  418.   110 GO TO 170
  419.  
  420.   120 IF (('unif').NE. (type)) GO TO 130
  421.       range = parin(2) - parin(1)
  422.       av = parin(1) + range/2.0
  423.       var = range**2/12.0
  424.       GO TO 170
  425.  
  426.   130 IF (('beta').NE. (type)) GO TO 140
  427.       av = parin(1)/ (parin(1)+parin(2))
  428.       var = (av*parin(2))/ ((parin(1)+parin(2))*
  429.      +      (parin(1)+parin(2)+1.0))
  430.       WRITE (*,*) ' A, B, AV, VAR ',parin(1),parin(2),av,var
  431.       GO TO 170
  432.  
  433.   140 IF (('bin').NE. (type)) GO TO 150
  434.       av = parin(1)*parin(2)
  435.       var = av* (1.0-parin(2))
  436.       GO TO 170
  437.  
  438.   150 IF (('pois').NE. (type)) GO TO 160
  439.       av = parin(1)
  440.       var = parin(1)
  441.       GO TO 170
  442.  
  443.   160 WRITE (*,*) 'Unimplemented type ',type
  444.       STOP 'Unimplemented type in TRSTAT'
  445.  
  446.   170 RETURN
  447.  
  448.       END
  449.